Measuring seismicity diversity and anomalies using point process models: case studies before and after the 2016 Kumamoto earthquakes in Kyushu, Japan

This paper reviews seismic activity in and around the Kumamoto region before and after the April 16, 2016, Kumamoto earthquake of M7.3 using statistical models such as stationary, two-stage, and non-stationary epidemic-type aftershock sequence (ETAS) models to examine seismicity anomalies. Our findings are summarized as follows. First, most of the earthquake clusters before April 2016 are explained by the stationary ETAS model, except for a few clusters of swarm activity, one of which was remotely induced by the 2011 Tohoku-Oki earthquake (M9). The non-stationary ETAS model describes changes in the rate of background seismicity of swarm activity. Second, we revealed seismic quiescence relative to the stationary ETAS model in the foreshock sequence from the M6.5 earthquake on April 14, 2016, and further in the aftershock activity of the 2000 M5.0 earthquake that occurred in the shallower extension of the M6.5 foreshock zone. Thirdly, the main-fault and two off-fault aftershock clusters of the M7.3 mainshock show different features, caused by static triggering effects of the mainshock and/or effects induced by fault weakening. Finally, the b-value increased stepwise over time during the entire period of foreshocks and aftershocks, the reason of which is explained. © The Author(s) 2017. This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. Open Access *Correspondence: [email protected] 1 The Institute of Statistical Mathematics, Tachikawa, Japan Full list of author information is available at the end of the article Background The 2016 Kumamoto earthquakes, including the mainshock (M7.3), were a series of shallow, strong earthquakes that occurred at and around 01:25 JST on April 16, 2016 (16:25 UTC on April 15), near Kumamoto City in Kyushu, Japan; the first foreshock (M6.5) occurred at 21:26 JST (12:26 UTC) on April 14, 2016. Several major earthquakes occurred inland in Kyushu at shallow depths (≤ 30 km) within a 100–200-km radius preceding the 2016 Kumamoto mainshock. In January and April 1975, two events of M6.5 and M6.3 (the Northwest Kagoshima-ken earthquakes) occurred successively at distances of 40 and 65 km, respectively, south of the 2016 Kumamoto mainshock. A shallow M7.0 earthquake occurred in March 2005 (the Fukuoka-ken Seiho-oki earthquake) off the northern coast of Kyushu and 110 km north of the Kumamoto event. More recently, in November 2015, a shallow M7.0 earthquake (the Satsuma-hanto Seiho-oki earthquake) occurred off the west coast of Kagoshima Prefecture, Kyushu, approximately 200 km southwest of the Kumamoto M7.3 earthquake. The aftershock activities of the former and latter M7.0 earthquakes were followed by clear quiescence in the sense that until the 2016 Kumamoto earthquakes, no event of scale M4 or larger occurred for more than 10 years and 3 months, respectively. In addition, the March 11, 2011, M9.0 TohokuOki mega-earthquake occurred approximately 1200 km northeast of the Kumamoto event. Preceding the Kumamoto earthquakes, shallow background seismicity in and around Kumamoto Prefecture had long been present inland in the Japanese Archipelago (Ogata 2017c). Kumamoto Prefecture lies at the western extension of the Japan Median Tectonic Line, where Page 2 of 22 Kumazawa et al. Earth, Planets and Space (2017) 69:169 a system of active faults forks in two directions at the Beppu–Haneyama Fault Zone. Specifically, in April 2016, a series of earthquakes ruptured segments of the Hinagu Fault and the Futagawa Fault to its north; these earthquakes occurred along the southern boundary of the Beppu–Shimabara graben, with epicenters moving from west to east over time. Therefore, this area is known for its high seismic, volcanic, and geothermal activity and the seismic activity in this region exhibits extremely diverse patterns. To quantitatively illustrate this variety in seismicity, we analyzed microseismicity in the Kumamoto region since 2010, before the occurrence of the first M6.5 foreshock. We then analyzed the foreshock sequence of the M6.5 event until the time of the M7.3 Kumamoto mainshock. We further examined how aftershock activity in the mainand off-fault zones differed regionally. Finally, we evaluated temporal and spatial changes in b-values in the sequence throughout the M6.5 foreshock sequence and the M7.3 aftershocks over a 2-week period. For these analyses, we use the epidemic-type aftershock sequence (ETAS; Ogata 1985, 1988) model and its extended versions (Kumazawa and Ogata 2013; Kumazawa et al. 2016a, b) for a variety of seismicity patterns in addition to a model for the magnitude frequency changes (Ogata 1989; Ogata et al. 1991). For each dataset, we compare the goodness-of-fit of the models for the best fit, as given in the following section. Models and methods The ETAS model and its extensions The dataset {(ti,Mi); S < ti < T } comprises the occurrence times of the completely detected earthquakes in a target interval [S, T] associated with their magnitudes. Therefore, the stationary ETAS model (Ogata 1985, 1988, 1989, 1992, 2006a, b; Utsu et al. 1995), is used as the baseline model, where the five constants μ, K0, α, c, and p are estimated from each dataset. The parameter μ denotes the background seismicity rate, and the other four parameters control the portion of the seismicity rate triggered by the preceding earthquakes. The uppercase character S is the starting time of the target analysis period to which the model is to be applied, and ti represents the occurrence time of the ith earthquake associated with a magnitude of Mi that is greater than the cutoff magnitude Mc. The history Ht indicates that both coordinates of such earthquakes occurred before time t. Note that when large earthquakes precede the time S, the history Ht should include the occurrence time and the magnitude records of those large earthquakes. (1) θ (t|Ht) = μ+ ∑ { i: S<ti<t} K0e α(Mi−Mz) / (t − ti + c), For example, such a history includes the mainshock and aftershocks before time S even if these are incompletely detected. The simplest alternative model for the case where the stationary ETAS is misfit to a dataset is a two-stage ETAS model that uses different parameter values before and after the change-point time. If we hypothesize a fixed change-point time, this model is compared to the single baseline ETAS model using the Akaike information criterion (AIC; Akaike 1973, 1992). However, if the changepoint time is estimated, the goodness-of-fit relative to the baseline ETAS model is compared using a modified version of the AIC (Ogata 1992; Kumazawa et al. 2010) because the maximum likelihood estimate (MLE) of the change-point does not satisfy ordinary large-sample theory (Ogata 1978). Using X-window-based interactive graphical software for data selection (TSEIS; Tsuruoka 1996) and for the mounted ETAS analyses (XETAS; Tsuruoka and Ogata 2015a, b; Ogata and Tsuruoka 2016), the seismicity in several focal regions was preliminarily explored to determine the fit or misfit of the ETAS model. For the detailed manual and related software source, see Ogata (2006a). If the extended cumulative curve of the first ETAS model for the earlier period overpredicts the empirical cumulative curve in the second period, we may suspect quiescence relative to the first period (relative quiescence). One of the possible physical reasons for relative quiescence is that the seismogenic source is covered by the stress shadow (Ogata 2006b, 2007, 2011a). In contrast, relative activation can often be triggered by a major external earthquake occurring outside the focal region at time τ, and the seismicity rate can then be expressed by where the last term represents the triggered rate change caused by the external earthquake occurring at time τ. In the present case, the M7.3 Kumamoto earthquake has a potentially static or dynamic triggering effect (cf., Hill and Prejean 2015) to off-fault activities in the Aso and Oita regions. Furthermore, the 2011 M9 Tohoku-Oki earthquake may have triggered microseismicity in the central Kyushu region prior to the Kumamoto earthquakes. The non-stationary ETAS model (Kumazawa and Ogata 2013, 2014) extends the ETAS model (Ogata 1985, 1988) to fit transient swarms, including induced seismicity, as discussed by Hainzl and Ogata (2005). The parameters for the background rate μ and the aftershock productivity (2) θ (t|Ht) = μ+ ∑ { i: S<ti<t} K0e α(Mi−Mz) / (t − ti + c) + Kτ / (t − τ + cτ )τ , Page 3 of 22 Kumazawa et al. Earth, Planets and Space (2017) 69:169 K0 of ETAS model (model (1)) are set to be time-dependent in the non-stationary ETAS model: where μ(t) and K0(t) are technically characterized by piecewise-linear segments connected at the occurrence times of the earthquakes. For the estimation, occurrence data are inverted into the optimal solutions of μ(t) and K0(t) under proper smoothness constraints with the assumption that the aseismic stress and changes in the triggering parameter are both stably estimated using an empirical Bayesian method. Kumazawa and Ogata (2013) justified this estimation procedure using simulated datasets. Furthermore, Kumazawa et al. (2016a, b) derived consistent μ(t) results for a number of volcanic swarms with another independent approach that employs hourly sampled volumetric stress changes as external records. The goodness-of-fits of the above extended models are compared based on the increments ΔAIC = AIC − AIC0 or ΔABIC = ABIC − AIC0, where AIC0 represents the AIC value of the baseline single ETAS model. Note that the baseline ETAS model is the same as that in the case of the reference ETAS model of the non-stationary model. Specifically, when the weights of the constraints are very large, μ(t) and K0(t) become the same as the constant values μ and K0 in model (1), respectively, as the corresponding reference parameter. The datasets used throughout the present paper are part of the Hypocenter Catalog compiled by the Japan Meteorological Agency (JMA 2017), and the legitimacy for the results of our analyses performed on those datasets is given in Discussion section. Magnitude frequency distributions Changes in b‐values in space and time Consider the Gutenberg–Richter formula (Gutenberg and Richter 1944): with the constants a, b, and β = b ln 10. Restricting the range of earthquake sizes such that M ≥ Mc, we can derive the probability density distribution f(M|β) = λ(M )/ Λ(Mc) = βe−β(M−Mc), where �(Mc) = ∞ Mc (M)dM is the expected total number of earthquakes with M ≥ Mc. The β value in Eq. (4) can depend on location or time such that β(z) is a function of location z = (x, y) or time z = t. As the maximum likelihood estimate of the b-value is given by the reciprocal of the magnitude average, local changes in the b-values have conventionally been obtained via various kernel methods or moving weighted (3) θ (t|Ht) = μ(t)+ ∑ { i: S<ti<t} K0(ti)e α(Mi−Mz) /